#在下载数据前的预设值
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GEOquery")
BiocManager::install('seandavi/GEOquery')
packageVersion("GEOquery")
BiocManager::install("Deducer")
library(ggplot2)
library(rJava)
library(JavaGD)
library(JGR)
library(car)
library(carData)
library(MASS)
library(Deducer)
library(Biobase)
library(BioGenerics)
library(GEOquery)
library(dplyr)
library(tidyverse)

#下载数据并处理
gset = getGEO('GSE38959', destdir=".", AnnotGPL = T, getGPL = T)#数据
exp<-exprs(gset[[1]])  
pdata<-pData(gset[[1]])
group<-c(rep("control",3),rep("hht",3))	
GPL<-fData(gset[[1]])
gpl<-GPL[,c(1,2,3)]
gpl$`Gene symbol`<-data.frame(sapply(gpl$`Gene symbol`,function(x)unlist(strsplit(x,"///"))[1]),stringsAsFactors=F)[,1]
gpl=as.matrix(gpl)
exp_symbol=cbind(gpl[,2],gpl[,3],exp)
exp_symbol<-na.omit(exp_symbol)
table(duplicated(exp_symbol[,1]))
exp_unique<-exp_symbol[!duplicated(exp_symbol[,1]),]
BCancer=exp_unique[,1:45]
n=dim(BCancer)[1]



#定义三个多重检验函数
#BH过程
BH_RD=function(ve,a,K){
nu=K
ve1=sort(ve)
c2=rep(NA,time=K)
for(i in 1:K){
c2[i]=(i*a)/K}
for(i in K:1){
if(ve1[i]>c2[i]){nu=nu-1}
else{break}
}      
ve2=c(nu,order(ve)[1:nu])
}
#e-BH过程
e_BH_RD=function(ve,a,K){
nu=K
ve1=sort(ve,decreasing=T)
c1=rep(NA,time=K)
for(i in 1:K){
c1[i]=K/(i*a)}
for(i in K:1){
if(ve1[i]<c1[i]){nu=nu-1}
else{break}}            #此时nu是拒绝的假设个数。
ve2=c(nu,order(ve,decreasing=T)[1:nu])
}
#e-up-down过程
e_up_down_RD=function(ve,a,K,p){
ve1=sort(ve,decreasing=T)
ve2=order(ve,decreasing=T)
c2=rep(NA,time=K)
for(i in 1:K){
c2[i]=K/(i*a)}
if(ve1[p]<c2[p]){
nu=p
for(i in p:1){
if(ve1[i]<c2[i]){nu=nu-1}
else{break}
}               
}
if(ve1[p]>=c2[p]){  
nu=p-1
for(i in p:K){
if(ve1[i]>=c2[i]){nu=nu+1}
else{break}
}         
}
return(c(nu,ve2[1:nu]))
}

#计算p值和e值并整合数据
library(dplyr)
E=array(data=0,dim=c(n,5))
colnames(E)=c("gene_name","gene_symbol","fold change","E-values(W)","p-values(W)")
for(i in 1:n){
E[i,5]=as.numeric(wilcox.test(as.numeric(BCancer[i,3:32]),as.numeric(BCancer[i,33:45]), exact=FALSE, alternative = "two.sided")[3])
E[i,1]=BCancer[i,1]
E[i,2]=BCancer[i,2]
E[i,3]=mean(as.numeric(BCancer[i,3:32]))/mean(as.numeric(BCancer[i,33:45]))
E[i,3]=log2(as.numeric(E[i,3]))
for(i in 1:n){
p=as.numeric(E[i,5])
if(p>=exp(-1)){E[i,4]=1}
else if(p<exp(-1))
{E[i,4]=-exp(-1)/(p*log(p))}
}

#表5.13数据代码(仅写出a=0.05的情况)
a=0.05
t1=BH_RD(as.numeric(E[,5]),0.05,n)
t1[1] #BH过程识别的基因数
pdata=E[t1[2:(t1[1]+1)],]
sum(abs(as.numeric(pdata[,3]))>1)/t1[1]#BH过程显著基因所占百分比

t2=e_BH_RD(as.numeric(E[,4]),0.05,n)
t2[1] #e-BH过程识别的基因数
edata=E[t2[2:(t2[1]+1)],]
sum(abs(as.numeric(edata[,3]))>1)/t2[1]#e-BH过程显著基因所占百分比

t3=e_up_down_RD(as.numeric(E[,4]),0.05,n,n)
t3[1] #e-down过程识别的基因数
edata=E[t3[2:(t3[1]+1)],]
sum(abs(as.numeric(edata[,3]))>1)/t3[1]#e-down过程显著基因所占百分比

t4=e_up_down_RD(as.numeric(E[,4]),0.05,n,(n+1)/2)
t4[1] #K/2-e-up-down过程识别的基因数
edata=E[t4[2:(t4[1]+1)],]
sum(abs(as.numeric(edata[,3]))>1)/t2[1]#K/2-e-up-down过程显著基因所占百分比




 

